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PREFACE 


This is Part IV of a seven part final report under 
Contract No. NAS 8-21143 between rhe George C. Marshall Space 
Flight Center and the University of Alabama, This report 
includes the results of a systems analysis of the transienrs 
occurring in an assumed model of the proposed flight experi- 
ment for boiling and vapor bubble studies. The technique of 
digital simulation is employed in the analysis. 

The basic equations of fluid mechanics and heat transfer 
are determined for each component making up the assumed 
flight experiment model. This includes the fluid flow circuit, 
the heat exchangers for boiling and condensation, the pump 
and the motor. These equations are solved simultaneously for 
a limited range of variables by digital simulation using a 
high speed computer. 

The results indicate that digital simulation can be a 
useful tool in determining the total system characteristics. 
This is expecially useful in determining start up or shut 
down transients as well as transients resulting from 
catastrophic events such as value failure. 


V 



CHAPTER I 


INTRODUCTION 

Early liquid propellant rockets experienced difficulty 
in maintaining a steady flow of fuel to the rocket: engines 
The tro\ible in these early rockets was due to the turbulent 
flow conditions existing at the interface betv^een the fuel 
feed lines and the propellant tanks. Fine mesh screens and 
thin porous blankets were found to be effective in reducing 
these undesirable effects. Recently the more efficient 
cryogenic liquid fuels are being used, a characneristiic of 
these cryogenic liquid fuels is a tendency no boil or 
vaporize in the propellant tanks during flight. This vapor 
collects on the above mentioned screens and obstructs the 
flow of the liquid fuel to the engines, A possible solution 
to this problem may be replacing the screens with coarsely 
packed porous beds. 

Therefore, in order to gain more fundamental knowledge 
in the area of two-phase flow through porous beds, the 
National Aeronautics and Space Administration in Huntsville, 
Alabama awarded the University of Alabama the contract to 
design a flight experiment to determine the behavior of 
two-phase vapor-liquid and/or gas-lic[uid flow through 
porous beds in a low gravity environment. This thesis is 
concerned specifically with that part of the experiment: in 
which the second phase (vapor) is obtained by heating the 
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liquid (water) to lioiling hy using electrical- thermal 
elements within the porous bed. 
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An idealized mathematical model of the proposed system 
is developed. The resulting system of differential and 
algebraic equations is digitally simulated using the 
International Business Machine (IBM) System/360 Continuous 
System Modeling Program (CSMP)« The objective of the 
simulation is to study and evaluate the modeling of the 
system? however, this thesis does not give an e:sihaustive 
study of the system but shows representative results 
obtainable using the simulation. 

In Chapter II the proposed system is divided into 
several subsystems and the equations of motion are developed 
for each subsystem. Chapter III begins with a discussion of 
some of the salient features of CSMP as they pertain to this 
thesis and concludes with the development of the simulation 
program including a flow-diagram. Within Chapter IV the 
discussion and display of the simulation results are given. 
In Chapter V, the final chapter in the thesis, concluding 
remarks with regard to the digital simulation and results 
obtained from the simulation are given. 



CHAPTER II 


DEVELOPMENT OF THE SYSTEM EQUATIONS 

In this development, the mathematical models derived 
are for an idealized system* An idealized system implies 
the ultimate in functional capability and is not limited 
by either physical or economic factors* The idealized 
model conveys the functional sense in which the system is 
to be used and is sufficiently complete, while avoiding 
attention to fine details* 

The development of the “equations of motion** for the 
fluid system is now given* The expression “equations of 
motion" is used for both differential and algebraic 
equations within the system. Assumptions and/or approx- 
imations made, are stated during each mathematical 
development. 

This system is a “mixed system" in that electrical, 
mechanical, fluid and thermodynamic quantities appear 
together and are, of course, interrelated, A diagram of 
the entire system is given in Figure 1, The system is 
discussed in this chapter under the following headings: 
Motor and Speed Control 
Positive-Displacement Pump 
Preheater and Control 
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FIGURE i. SCHEMATIC DIAGRAM OF TPIE FLUID SYSTEM. 
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Vaporization Chamber and Condensation 
Channel and Porous Bed 
Accumulator and Vaporization Control 
Pressure Equations 
Motor and Speed Control 

A linear DC-motor using field voltage as a means to 
vary speed is used. It is assumed the motor is supplied 
with a constant armature current. 

This means for regulaning the speed of the DC-motor 
is realized through the use of state-variable feedback. 

Conventional methods of synthesizing a linear control 
system usually call for some type compens anion network to 
be designed to achieve the performance characteristic 
needed to perform the syscem’s designated cask. Through 
the use of state- variable feedback, the task of synthesizing 
a compensating network is reduced to that of determining the 
appropriate feedback coefficients. This is done using 
simple algebra and block diagram manipulations. The 
designer selects the time response and then Laplace 
transforms it into a rational fraction of the complex 
variable s. The system is depicted in standard block 
diagram form where all the state variables are fed back 
through the appropriate transducers to the inpuc of the 
system. In following the above procedure, a system can 
be made to conform to given time domain specifications 
as long as those specifications can be translated into 
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the appropriate mathematics.^ 

The state- variable synthesis technique is now 


demonstrated ; 


,o 

l5if + Rif = Vf 

(1) 

JW +BW=T-Tl-Td 

(2) 

Vf s= Ka ERROR 

(3) 

T = K^if 

(4) 

ERROR REF - kiW - k2if 

(5) 


Equations (1) and (2) have been time scaled so that the 
unit of time employed in this development is minutes. 

t = time in minutes 

if = the field current of motor (amp) 
o 

if s= dif/dt (amp/min) 

iJ « the field inductance (henrys)/60 

R = the field resistance (ohms) 

o ’ 

W « the derivative of angular velocity with 

respect to time 

J = the polar moment of enertia of the motor plus 
the piimp (lbf-in2/386.4) 

Ibf = pounds force 

B “ viscous friction of the motor plus the pump 
( Ibf-in-min) 

^Superscripts denote references given at the end of 
this thesis. 
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T e the electromagnetic torque of the motor 
(lhf«in) 

Tl ® the torque needed to overcome the fluids 
resistance to flow 

Td b any arbitrary disturbance torque applied to 
piotor shaft C3.bf-"in| 

Vf = the applied field control voltage (volts) 

Ka ® the DG~amplifiar gain 

Km - the electromagnetic torque constant, of the 
DC~motor (Ibf-ln/amp) 

REP = the applied reference DG-=voltage to the motor 
speed control (volts) 

ki = feedback coefficient (volts<-min/rad) 

k 2 = feedback coefficient (ohms) 

A block diagram representation of the Laplace 
transformed electrical and mechanical Equations 1 through 5 
is given in Figure 2o 

The motor and the mechanical load forms a seeond~order 
system® A second-order system which is sl,lghtly under- 
damped is assumed to be adequate® Second-order systems 
have been thoroughly studied and there exist normalised 
plots^ for different undamped natural frequencies (Wji) and 
damping factors (Sn)° Therefore) the desired time response 
was easily transformed into a rational fraction of the 
complex variable s and is given byg 

6 ( s ) s= Wn^/^ ^ ^ ® ) 

where the coefficients and completely determine the 
response of Equation 6 to a unit-step input® 



FIGURE 2„ BLOCK DIAGRAM OF THE MOTOR SPEED CONTROL SYSTEM 
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In order -to consider -the system shown in Figure 2 to 
be a single-input, single-output system, the torques Tl 
and Td are set to equal to zero. 

The block diagram of the speed control system in 
Figure 2 is now manipulated into a convenient form so the 
feedback coefficients may be determined „ This is accom- 
plished in two steps* See Figures 3(a) and (b) * 

Equating the results of figure 3(b) to Equation 6 
gives the results in Equation 7* The feedback coefficients 
1^1/ ^2 and the gain Ka may be determined directly by 
equating like coefficients of the left and right sides of 
Equation 7 * 



t 2SiiW ns 



(7) 


Ka Km/1? J 

— I II etj i o tii M iia L ii i Ci ^ Ii i ■ I I- II r m ii i i r^ »i i rtf i f i — 

s^ «}• (R/L t B/J t k2KaJ)s + (K^Kj^kx + RB/LJ t k2KaB) 


where s 

Wn = the known undamped nattjral frequency of a 
secopd-order system (rad/min) 

Zxi - the kno\«i daartping factor of a second-order system 

After kl, k2 and ka are determined, the motor speed 
control system as shown ip. Figure 2 in mathematically 
equivalent to Equation 6 from the input, output viewpoint. 
Posit! ve-Displacemen'c Pump 

Before beginning the development of the positive- 
displacement pump, it is well to note a couple of general 



REF 




FIGURE 3, BLOCK DIAGRAM REDUCTION OF SPEED-CONTROL SYSTEM 


o 
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characteristics of this type of device. 

Probably the most distinguishing feature of this 
device is that the quantity of liquid displaced per 
revolution of the shaft is approximately a constant. 
Another important characteristic is that the pressure 
developed by the pump is a function of the fluid circuit. 

This mechanical-fluid device is analogous to the 
constant current generator from electrical circuit 
theory if pressure is the analog of voltage and fluid 
flow is the analog of current. ' 

An ideal positive-displacement pump is defined as 

having no power losses due to friction and leakages and, 

consequently, has an efficiency of 100%. In practice, 

80% efficiency is readily obtained and, therefore, as a 

first approximation, the system analyst will often use 

. 3 

the ideal pump as the model for the actual devxce. 

Consider now the equations for the ideal positive- 
displacement pump: 

Mechanical Power In = Fluid Power Out 
W = A P Q 

1j 

The torque, T^ (Ibf-in) , is required to overcome the 
back pressure of the fluid circuit in order to produce a 
given flow rate. The back pressure existing in a fluid 
circuit is a function of the physical parameters of the 


circuit. 
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The angular speed of the pump shaft is W. Since 

the motor is directly coupled to the pump, the motor speed 

and the pump speed are the same (rad/min) , The differential 

2 

pressure drop across the pump is AP (Ibf/in ) . 

The volumetric flow rate (Q) is normally expressed 
linearly in terms of the displacement constant (Dj^) and 
the angular speed of the pump shaft. 

Q = W (9) 

4 

Factors affecting the accuracy of Equation 9 are: 

1. Pump speed 

2. Pressure at discharge and pressure on suction 
side 

3. Viscosity of fluid 

4. Amount of entrained air or non- condensable 
gases in fluid 

of the above mentioned factors , the two most important 

for this system are numbers (2) and (4) « 

The pump used is assumed capable of delivering five 

gal/min into a back-pressure of at least an order of 

magnitude higher than is expected to occur during the 

experiment. The effects of large pressure differentials 

across the pump are to increase the amount of slip flow 

which decreases the efficiency of the device. Slip flow 

is defined as the flow in the reverse direction through 

, 5 

a pump due to large pressure differentials. 

If two-phase flow were to occur through the pump, its 
exact effect would be difficult to predict, but it may be 
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concluded that the volumetric constant of the pump would 
effectively be reduced. This conclusion is drawn from an 
intuitive concept that since displacement volume would be 
shared by fluid and vapor, there would naturally be less 
fluid delivered per revolution. Therefore, to deliver the 
same flow the pump speed would have to be increased. 

As implied from the above discussion, the entrainment 
of vapor through the pump is a highly undesirable phenom- 
enon. It is assumed that this phenomenon does not occur. 
The validity of this assumption is checked with the 
simulation. 

Preheater and Control 

The first law of thermodynamics is sufficient to 
describe the process of heating flowing water in a tube, 
a typical cross-section of the preheater which is shown 
in Figure 4 . 

! rate of change of energy > _ 

Within the control volume/ 

^ rate of heat flow-in + rate that heat\ 
vis added by the heaters J 

/rate of heat flow-out + rate of heatN 
V transferred by convection / 

It is assumed the tubing and water are at the same 
temperature. Also transfer of energy from one section of 
water to another by conduction is neglected since the flow 
rates are much larger than the heat conduction rates . • Heat 
conduction would be significant if the system was operated 
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statically (no flow) , with the heaters on since conduction 

would be the only thermal coupling between sections. But 

by the very nature of the experiment, the system will 

always operate dynamically. 

The equations of motion of the preheater are given 

as follows : 
o 


P V C T. . 
w w w wx 

o 


P C Q 
w w 


(T 


w (i-1) 


’’wi> + 


P V e T = UA (T . 
a a a a wx 


T^) + - T^) 


o 


( 11 ) 


PjV^CjTj = (T, - T^) 


1 " 1 ' 


U2*2<’'o - 


(12) 


where : 
o 

T = dT/dt 

3 

P^ = the density of water (Ibm/in ) 

3 

V = the volume of water in the ith section (in ) 
w 

C = the specific heat of water (BTU/lbm°F) 
w 


Q = the volumetric flow-rate of water (in'^/min) 

U = the heat transfer coefficient between the 
tubing and air (BTU/in^ min °F) 

A == the heat transfer area for the tubing and 
axr (xn ) 

Ql = the heater input (BTU/min) 

3 

P = the density of air (Ibm/in ) 

Ci 

3 

V = the volume of air in the ith section (in ) 

a 

= the specific heat of air (BTU/lbm °F) 

U, = the heat transfer coefficient between air 

^ 2 o 

and insulation material (BTU/in min ^P) 
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A- = the heat transfer area between the insulation 
and air (in^) 

3 

Pj =,the density of the insulation (Ibm/in ) 

V 3 

I = the volume of the ith section of insularion (in ) 

Cj. = the specific heat of insularion (BTU/lbm °F) 

= ambient temperature (°F) 

The preheater control is an "on-off" or "bang-bang" 
type. In this particular case, the heat-input (Ql) in 
the preheater section is either completely on or completely 
off. Implementation of the "on-off" control is through the 
control law. The control law in this case consists of 
several logical statements which are: 

If T^^ ^ (Reference Temperature - DX) then Ql - 0.0 

If T < 

wi — (Reference Temperature - DX-BAN) then Ql = Qmax 

If T . < (Reference Temperature - DX) , and Tv^i is 
positive in sign, then Ql = Qmax 

If T , < (Reference Temperature - DX) but > 

(Reference Temperature - DX-BAN) , and 

T^^ is megative in sign, then Ql =• 0.0 

An arbitrary increment of temperature, DX, separates 
the upper limit of the BAN and the Reference Temperature 
(saturation temperature as determined by the absolute 
pressure in the ith section of the tubing) . To account 
for inaccuracies in the control system, the increment, DX, 
was included to ensure the temperature T^^ would not reach 
the reference temperature. 
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The arhitrary increment of temperature iis BAN in 
which the heat-input (Ql) can bo off« Whefirtar the 
heat-input (Ql) is off within the BAN incremsvft is 
determined by the sign of the derivative of tiie water 
temperature (dT^^i/dt =» ^wi^ within the ith section of 
tubing# The width of this temperature increment is set 
by the requirements on the amount of vaporization 
needed to successfully accomplish the flight experiment, 

A time response like that shown in Figure 5 is 
representative for the preheater control system, 

Xf vaporization occurs in the preheater section. 

Equation 10 must be replaced by another equation in 
order to account for this change of phase. The derivation 
and explanation of this phenomenon is given in a later section. 
Vaporization Chamber and Condensation 

The saturation temperature {boiling point of water) 
is a function of both pressure and temperature, Forma'fcion 
of vapor bubbles is considered negligible at temperatures 
below saturation temperature. At saturation tempez-ature; 
the vapor pressure tends to exceed the local absolute 
pressure which results in the development of vapor bubbles. 
Enough vapor bubbles are formed to prevent the temperature 
of the water to rise above the boiling point. Therefore^ 
the temperature of the boiling water remains constant 
until either all the liquid is vaporized or there is a 
change in the absolute pressure. 




Q1 (BTU/MIN) (Twi°F) 
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FIGURE 5, TIME RESPOHSE OF THE PREHEATER CONTROJ^ 
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In this presentation all the boiling is asstuned to 
occur in the “vaporization chaiiiber” even though localized 
vaporization due to extreme pressure gradients may take 
place , 

hn energy balanqe is written for the cross-section of 
the vaporization chamber in Figure 6, 


B rate of energy stored \ 

\within the liquid / 

f rate of heat flow-in rate 
Vthat heat is added by the heaters 

^rate of heat flow-out -i- vapor boil-up rate 

d/dt(P„V;/:„Twc) “ •{- 


PwCwQ(^W”'^wc) UA(To ~ ^wc) 


(13) 


Upon comparing Equation 10 and 13 j two differences are 
notedo Firsts in Equation 13, the derivative of the entire 
terra on the left side is taken since the volume of water is 
changing with vaporization. Second, the term 'X v appears 
on the right side of Equation 13, This term accounts for 
the loss of heat due to vaporization where v is called 
vapor boil-up rate and X is the latent heat of the vapor. 

Performing the differentiation of Equation 13 gives 


PviC<MTv7cVw J*w^w^w^wc ~ Q2 

•f’ PvjCvjOCTi^jf — T^q) + UA(Tq — T^(^) — ^v, 

o 

Substituting Pw^w ® into Equation 14 gives 


( 14 ) 
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V « PvjCwQ/(7\ —Cw^wc) *" ^wc) 

UA3/(A —C^T'^q) (Tq " ^V7C^ ** (15) 

^wCw^wTwg/(A —C^T v^c' Q2/(A “G^Ti^c) « 

Now Equation 15 gives the rate of vaporization in terms of 
fluid temperature and the heat input Q2o The method of 
controlling the vaporization will he given in Section VI « 
Condensation is the negative of vaporization and is 
dependent upon the heat losses in those sections in which 
vaporization occurs* Recalling the assumption that all the 
vaporization occurs locally within the vaporization chamher^ 
hut condensation will occur in sections of the system which 
contains vapor huhhles* It is assumed that vapor hubbies 
cannot exist in those sections in which the surrounding 
fluid temperature is below saturation. These assumptions 
are discussed in a later section with regard to a modified 
definition of thermodynamic quality. The equauions describ- 
ing condensation within the system are of the same form as 
the vaporization expression, Equation 15, There is a 
condensation equation for each section of the system in 
which vapor bubbles exist, A typical expression is 


CONDEN = PwCwQ/(^ -CwTwc) (Twc - ) + 

UA3(T^(j — *I’o)/( ^ “^w^w) 


( 16 ) 


where the terms in Equation 16 are similar uo those in 
Equation 15, 
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Channel and Porous Bed 

The porous bed consis-bs of small glass spheres packed 
in the channel o How well -the porous media is packed can 
be expressed by iihe so called ’’porosiby of t.he bede” 

Figure 7 shows a section of the bed displaying the porous 
media* Porosiry { 0 ) is defined as the ratio of the pore 
volxime to the total volume 

The overall objective of this experiment is to 
determine the behavior of two-phase vapor-liquid flow 
through porous beds in low gravity environment* Therefore, 
to the fluid dynamist and the chermodynamist., the phenomenon 
occurring within the channel and porous bed is the focal 
point of their work* But to the system analyst, the channel 
and porous bed form but a lirik of an entire system* In keep- 
ing with the system® s point of view, the equations of motion 
for this section are derived by simplifying the problem as 
much as possible* 

Equation 17, which had been developed by another 
investigator, is used as a starting place in this develop- 
ment 


dP/dX =C<u Vi/d2 + /pu Vi/d2^ Re 


(17) 


Equation 17 allows one to calculate the pressure at any 
point along the length of the porous bed as a function of 
the experimental parameters 0< ^ p; and the Reynolds Number 
(Re) for the flow* Reynolds Nx^mber in the chanrAei is 
defined to be the ratio of the product of fluid density. 



ChanHiel 


Pore space 



FIGURE 7o CROSS-SEGTiOmL VIEW OF CIIAENEL SHOWING POROUS 
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the diameter of one sphere, and the average fluid velocity 
to the fluid viscosity. The other parameters in Equation 17 
are d, the diameter of a sphere, the velocity of the water^ 
and u, the viscosity of the water. The temperature and 
pressure are no longer independent properties of water once 
saturation is reached. So, if the pressure is known; the 
temperature of the saturated water is automatically 
determined. With the dependence of temperature upon pressure 
along with results of Equation 17; the temperature variation 
of saturated water in the porous bed can be determined. 

During times when the water within the porous bed is not 
saturated, then an equation of the form of Equation 15 with- 
out the terms 7 \ v and Q2 is applicable for the unsaturated 
water. 

In order to make the analysis of the channel and porous 
bed complete, expressions for the temperature of the porous 
media and channel as they are coupled to the water 
temperature must be developed. 

The porous media is assumed to be equally distributed 
in a cross-section and along the length of the bed. From 
the knowledge of the volume of the bed; porosity of the bed, 
density of porous media, and diameter of one sphere, the 
total surface area and the total mass of the porous media can 
be calculated. If the bed is divided into n sections, then 
the surface area of the porous media in the nth section is 
the total surface area divided by n, the number of sections. 
Similarly, the mass of the porous media in the nth section 
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is the total mass of the porous media divided hy the 
numher of sections. The surface area of the porous media 
is needed in order to calculate the heat transfer from the 
water to the porous media or vice versa. In order to 
calculate the heat capacity of the porous media^ the mass 
of the porous media is required. 

The following procedure is used to obtain the area 
and mass of a sections 


Vs 

= (l-^)Vc 

(18) 

Ms 

II 

(19) 

Bn 

= (l-i^)Vc/4/3lVr| 

(20) 

Ba 

= 3(l-;5)Vc/rb 

(21) 


where s 

Vg = the volume of porous media in channel (In^) 

Vc * "the volume of channel (in^) 

Mg « mass of porous media in channel (Ibra) 

3 

^*pm = ■the density of the porous media (Ibm/in ) 

Bn “ "tli® number of spheres in channel 
rb - "the radius of a sphere (in) 

2 

Bg = the total surface area of the glass sphere (in ) 


Now the mass of the porous media in section n is 


Ms(n) “ BpjijVg/n. 


( 22 ) 
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The surface area of the porous media in section n is 

Ba(n) = Vc/rb n» 


( 23 ) 


Now the energy equations can he written for the nth section 
of the porous rtiedia^ the water in the nth section of the 
porous hed, and the nth section of the channalo See 
Figure 8» 

The energy balance for the nth section of the channel 

can be expressed as follows? 

/ rate of accumulation of energy j ^ 

\ within section n of channel / 

(^energy transferred from the water to the channel - 
energy transferred from the channel to the air ^ 
Fc^c^cl'cCn) ~ ^3^3 (“^wCn) “ ^c(n)) 


( 24 ) 


U4A3(To “ ‘2^c(n)) 


where? 

Pc - the density of channel material (Ibm/in^) 

Cc = the specific heat of the channel material 
(BTU/lbm~°F) 


Vc - the shell volume of the nth section of 
the channel (in'^) 

U3 = the heat transfer coefficient between uhe 
water and channel (BTU/in^-min-^F) 

A3 = the heat transfer surface area (in^) 

U4 = the heat transfer coefficient between the 
channel and air (BTU/in 2 ~min-®F) 

The energy balance for the nth section of the porous 


media is 
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I rate qf accumulation of energy 
(^within the nth section of porous media 

/ energy transferred from the water j 
\to the porous media in section n / 


PgVgCgTg = ^5^4^^w(n) " *^g(n)^ 


(25) 


where s 

Pg s= the density of the porous media (Ihm/in^) 

Vg * the volume of the porous media in the nth 
section (in*^) 

Cg “ the specific heat of the porous media (BTU/lhra-°F) 

U 5 = the heat transfer coefficient between the porous 
media and water (BTU/in^-min-^F) 

A 4 = the heat trart^fer surface area for water 
and the porous media (in^) 

Equation 26 describes the temperatures behavior of the 
water at the nuh section of the channel. The development 
of Equation 26 is similar to Equation 10 and will not be 
given, 

/ rate of accumulation of energy in the \ _ 

\nth section of water in the porous bed/ 

/ heat-in due to flow + heat transfer from the ) _ 

porous media ^ heat transferred from channel/ ” 

( heat-out due to flow j 

^wf (n)^wf (n)Cwf (n)^wf (n) ~ 


( 26 ) 


Pwf (n)Cwf (n)Q(T'wf (n-1) - ) + 

U5A4(Tg(n) - Tvrf(n) ) + U3As(Tc(n) “ Twf(n) )• 
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The symbols in Equation 26 are similar to those which 
have been previously definedo 

A summary of the results obtained in this section is 
now given. The channel was subdivided into several 
sections. The porous media was assumed to be evenly 
distributed throughout the channel. Then the total 
surface area and total mass of the porous media was 
determined. The surface area of any section is equal to 
the total surface area divided by the mamber sections, 
n. The mass of the porous media in any section is 
determined similarly. The time-dependent equations for 
temperatures of the channel, of porous media, and of the 
water were written. Thus the mathematical model for the 
process within this section was obtained. 

Equations 24, 25, and 26 are valid as long as uhe 
water temperature is below saturation. After saouration 
has been reached. Equation 26 is replaced by Equation 27 
with Equations 24 and 25 remaining the same. The appropri- 
ate equation was employed according to logical rules 
programmed into the simulation and is discussed in Chapter 
III. 


Twf(n) F(Pn) (27) 

The function, F(Pjj), expresses the relation 
between the temperature, Tvjf(n), in the nth section of the 
channel and the absolute pressure, in the nth section 
of the channel when saturation has occurred. The function 
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F(Pn) is constructed from data taken from thermodynamiq 
tables of saturation pressure and temperature* 

Accumulator and Vaporization Control 

The accumulator and vaporization control are presented 
together because in this system they are integrally related. 

Accumulators are widely used in fluid systems to smooth 
transients such as line shock; valve -surges, and pump flow- 
ripple, Any time a positive-displacement pump is used, as 
in this system, there will be some flow ripple. The amoxont 
of ripple depends upon the particular pump,® 

A simple example of an accumulator is an open tank with 
a piston or weight sitting on top of the fluid. See Figure 9. 
The piston is free to move up and down depending upon the 
level of water in the tank, 

f rate of accumulation) f 

water in tank J = { - Flow-out ; 

V = Qln - Qout (28) 

The pressures at the entrance and exit of the accumu- 
lator as given in Figure 9 ares 


P4 = Wp/Ap 

(29) 

P5 = Wp/Ap +YwT^ 

(30) 


where 8 

Wp s= the weight of the piston (Ibf) 
Ap « the area of the piston (in^) 
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h* *= the height of the piston above the input 
flow (in) 

h = the height of the piston above the output 
flow (in) 

K the specific weight of water (Ibf/in*^) 


The volxme of the accumulator as shown in Figure 9 is 
V = Ap h. Assuming that Ap and Wp remain constant; then 
Equation 31 follows by differentiating Equation 30 with 
respect to time. 


o 

PS 



= Vw/Ap (Qin-* Qout) 


(31) 


Equation 31 can be integrated to give the absolute 
pressure of the accumulator as a function uime<. If 
vaporization occurs in the system, then Equation 31 is 
modified to account for this effect « 

In this experiment; all foreign gases will be purged 
and the fluid circuit completely filled with water* This 
is done to ensure the bubbles formed are composed entirely 
of water vapor* When vaporization occurs and the vaporiza- 
tion minus the condensation is a positive quantity, the 
resultant effect on the system is an increase in uhe volume 
of the accumulator* It is assumed that water is essentially 
incompressible* The system is assumed to have an accumulator 
which has a component equation similar to Equation 31 and is 
given. by# 
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Pac = Kac Vac (32) 

Pac ” ^ac Vac (33) 

Vac == 1/Pwc (Qin - Qout) •»* SPV (VAPOR~GOE®EN) (34) 

where: 

K^c “ pressure~voluine constant (Ihf-min/in^) 

Pwc ” density of the water in the accumulator 
( lbin/in3 ) 

SPV = the specific volume of saturated vapor 

(in3/lbm) 


Therefore, ultimately, a net increase in the accumulator 
volume will manifest itself as a change in the absolute 
pressure around the fluid circuit® 

The vaporization control is motivated by two factors® 
First, results from initial simulation runs showed that the 
“On-Off" preheater control has marked effects on the 
vaporization during the “off" period of the preheater control® 
That is; the vaporization level will follow the variations 
in temperature of the water in the preheater section® Ideally 
a vaporization control will tend to smooth out the above 
mentioned variations in vaporization level. 

The second motivating factor is the need to continuously 
measure either quantitatively or qualitatively the thermo- 
dynamic quality of the two-phase vapor liquid® Quality is 
defined as^ 
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A modified definition of quality is also used in this 
developmento The denominator of Equation 35 is the total 
mass of liquid in the system where the denominator of the 
modified quality is the mass of the liquid herween the 
vaporization chamher and the accumulator* The numerator 
tern remains the same* The motivation for the alteyhate 
definition is that the experiment is to he designed so 
that vapor exists only between the vaporization chamber 
and the accumulator. When the term QUALl is used later 
in this thesis, it will denote the quality as defin-sd by 
Equation 35, while the term QDAL2 denotes the modified 
quality as defined above. 

The coupling between the vaporization control and 
quality measurement is not obvious. The mass of vapor in 
the system at any time is given by 

Mass of Vapor = I (Vaporization - condensation) dt.(36) 

Equation 36 can be integrated directly since vaporization is 
given in closed form by Equation 15 and condensation is given 
in closed form by Equation 16, Since in practice, vaporiza- 
tion and condensation are not readily measurable? therefore, 
another method of obtaining the mass of the vapor is needed. 
Equation 34 depends upon the flow-in, the flow-out, and the 
vaporization minus the condensation. The flow-in, flow-out, 
and volume, Vac» can be easily measured. The absolute 
pressure of the accumulator, Pac/ can be measured accurately 
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and it thus seems plausible to calibrate a system to give 
mass of vapor in terms of absolute accumulator pressure 
or volume* So any control scheme that used one of 

the control parameters would in a sense be -a thermodynamic 
quality control system. In the simulation it is assumed 
the liquid flow-in equals the liquid flow-out of the 
accumulator, since the overall vapor content is small as 
compared to the liquid. 

The proposed vaporization control is a pressure 
control system with the vaporization being controlled 
indirectly. It is noted there is no active way to control 
condensation other than manipulate the physical parameters 
of the system, A block diagram of the vaporization control 
system is given in Figure 10, The values for the feedback 
coefficients and A2 and the gain K^v can be obtained 
from the simulation. 

Pressure Equations for Fluid Circuit 

In this development it is assumed that pressure losses 
or pressure drops can occur in the tubing connecting the 
components themselves. Losses due to fittings; curature of 
the tubing, and instrumentation are neglecued as first 
approximation in this anaylsis. 

Tube pressure drop equations have been derived 
experimentally and are functions of the length (L) and 
diameter (D) of the tubing, viscosity (v) of the liquid 
and Reynolds Number (Re) of the flow. Equation 37 describes 
the pressure drop (Ap-t) in the tubing and it also accounts 




FIGURE 
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for entrance and exit losses. The entrance and exit losses 

are those incurred when the water is flowing from a smaller 

10 

diameter tube into a larger diameter tube or vice versa, 

A?t “ 128 u L (l-i-.0434 D/L Re) (37) 

Referring to Figure 9^ and Equations 29 and 30, the 
pressure drop, R 45 , from entrance to exit of the accumu- 
lator is given by 

P45 “ ^4-^5 ( 38 ) 

P 45 = X w<^’ ” (39) 

where (h® - h) is the difference in height from entrance to 
exit. It is assumed the height difference in the accumulator 
is negligible. 

Most of the equations presented thus far have been 
differential equations with time the independent variable, 
but the pressure equations, with one exception, are strictly 
steady-state. The exception being that the absolute pressures 
around the fluid circuit are affected by the pressure Pac(t) 
which is a function of time. 

Equation 17 gives the pressure drop across the channel. 
Note that Equation 17 and 37 are of the same basic mathema- 
tical form. 

The fluid circuit is idealized into the schematic as 
shown in Figure 11, From this schematic, the pressure drops 
and the absolute pressures, as referenced to the accumulator 



Porous 

Pi Tubing bed 




FIGURE 11 „ PRESSURE SCHEMATIC OP FLUID CIRCUIT „ 


pressure, may be calculated for various parts of the fluid 
circuit. 
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PlO 

St 

Pl2 + P23 - P3ac Paco 

(40) 

P2ac 


P23 + P3ac 

(41) 

Pi 

=; 

Pl2 ’*• P23 + PSac Pac 

(42) 

P2 


P23 PSac Pac 

(43) 

P3 


PSac Pac 

(44) 

Po 


~Poac Pac 

(45) 

where s 

PlO 


t-he total pressure drop around the fluid 
circuit (Ibf/in^) 


P2ac 

s: 

the pressure diffe.reritial between locations 
2 and ac .In fluid circuit 

Pi 

2= 

the absolute pressure at location 1 in the 
fluid circuit <lbf/in^) 


P2 

22 

the absolute pressure at location 2, in the 
fluid circuit (Ibf/in^) 


P3 

2? 

the absolute pressure at location 3 in the 
fluid circuit (ibf/in^) 


Po 

12 

the absolute pressure au location '0 in the 
fluid circuit (Ibf/in^) 
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CHAPTER III 

DEVELOPMENT OP THE SIMULATION PROGRAM 

Originally, the selection of digital simulation was a 
personal preference of the author but since the system as 
described in the previous chapter requires the use of 51 
integrators and 58 function generators, digital simulation 
was necessary because of the limited amount of analog 
equipment available. An IBM developed program called 
“Continuous System Modeling Program" (CSMP) which runs on 
the IBM System 360/ Model-50 computer was used for the 
entire simulation. 

No attempr. to compare digital and analog simulation 
or to compare other digital simulation techniques with 
CSMP will be made since the author’s practical experience 
in the general area of machine computation is not broad 
enough to make such critical comments. Rather, a sximmary 
of the essentials on how to run a typical CSMP program 
will be given. For those desiring a more functional 
knowledge of CSMP, the IBM literature on this program is 
a good source 

Before beginning the discussion of the development 
of the simulation program of this thesis, a brief descrip- 
tion of CSMP will be given. The purpose of this 



description is not detail instruction in the use of CSMP 
but to make familiar those features of CSMP as related to 
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this thesis, 

CSMP is essentially a program for solving a system 
of ordinary differential equations, A system of equations 
of the most general class or ordinary differential equations ^ 
i,e,, non«linear with time varying coefficients can be 
programmed. Therefore^ it is apparent that CS&IP is capable 
of solving a large class of engineering problems. CSMP 
will have particular significance to those people concerned 
with design and analysis of dynamic systems. 

It is assumed in the following discussion of the 
essentials of programming a problem for CSr-iP that the 
mathematical modeling of the system has been done and a set 
of differential and algebraic equations of uhe fom of 
Equations 46 and 47 are to ba prograromed. 

Oik) (46) 

0 = a(t)lj XiJ + Klj (47) 

1 ^ , e o 0 n 

1^ 2ff o ft , D m 

1^ 2^ o • o 4 p 

= independent variable 
= dependent variable 
s= forcing functions 


i « 
j - 

K - 

wheres 

t 

Xi j 
Uik 
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» time-varying or constant coefficients 

Kij = constants or parameters which could stand 
alone in a mathematical expression 

A general outline of a CSMP program is given in 
Figure 12# The remainder of this discussion will refer to 
elements within this figure » Figure 12 contains six 
sections of which all or some of the six sections will 
appear in every CSMP program# A CSMP program will always 
contain a set of CSMP control cards and the cards SND, STOP, 
AND ENDJOB, The other sections may or may noc he included 
in a particular program depending upon the particular 
problem at hand# 

Frequently, a program will contain parameters which are 
defined in terms of more basic system parameters and con- 
stants# Rather than perform these calculations before the 
program is run, they may be calculated within the INITIAL 
section# Calculations performed within the INITIAL section 
will be performed only once at the beginning of the simula- 
tion# In fact any calculation which needs to be performed 
only once at the beginning of the simulation can be done in 
the INITIAL section# 

Sometimes certain calculations are needed only at the 
end of the simulation# This may be achieved within the 
TERMINAL section of the CSMP program# An example of the 
use of the TERMINAL option would be parameter optimization 
in which it may be desirable to return to the main simula- 
tion program and rerun the simulation with the updated 




FIGURE 12 o BASIC STRUCTURE OF A CSMP PROGRAM 
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parameter values. When the RE-RUN opuion is used, the 
program x-zill rer.urn to the INITIAL section (if o’ue is 
included) where the new parameter value, as calculated xn 
the terminal section, will update the simulation <. The 
entire simulation will be rerun using the updated parameter 
value , 

A DYNAMIC section will probably be included in every 
simulation as the system differential equations are inte- 
grated within this section. This section is the most 
flexible withiii the program and can accommodate a wide 
variety of dynarai,c systems, i«e«, a wide range of mathe- 
matical formulations. 

An important feature of the CSMP program is the SORT 
and NOSORT options. If 5 section of coding contains no 
Fortran branching or logical tests then r.ne SORT option 
can be employed, Withlrx a SORT section, tiie program 
automatically places the statements In the bast possible 
sec[uence to ensure an accurate solution. Incorrect sequenc- 
ing of the structure sx^atements introduces a phase lag that 
could seriously affect the stability and accuracy of the 
solution. Within a NOSORT section, the burden of correctly 
sequencing the st'racture statements lies xvith the programmer. 
The advantage of uti listing the NOSORT option is that the full 
power of standard Fortran prograirtming is available. 

The discussion nov7 returns to the development of the 
simulation program for the system equations developed in 
the previous chapter. 
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The simulation program was developed in six steps „ 

These steps given in chronological order ares 

1, Separate or decouple the system equations into 
three independent sections s Electrical, 

Pressure and Thermodynamic* 

2* Assume constant values for all coefficients of 
the differential equations* 

3* Make appropriate runs on each section to check 
them out* 

4* Consolidate the three sections and make appropriate 
runs to check out system* 

5* Relax the constraint of step two where applicable 
and make appropriate runs* 

6. With steps one through five completed, make 
simulation runs with all systems functioning* 

The above six steps are not intended no imply any 
generality of method but to convey the method of approach 
taken in this problem* In general, if the system to be 
simulated is large, it is recommended that it be simplified 
and broken down into smaller subsysfiems as much as possible* 
This procedure is invaluable in the initial check-out and 
debugging phase of the program* 

The only equation coding problems occurred in what 
could be called “equation control*" Equation control is 
defined in this discussion to be the sv/itching-in and 
switching-out of the appropriate set of equations during 


the simulation* 
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The fluid-thermodynamic portions of the system operate 
either in the unsaturated state or the saturated state* In 
Chapter II, the fluid-thermodynamic equations of motion 
were developed. Though the exact details of these equations 
will vary depending upon the particular section of the system 
under investigation, they will he of the general form, 

PVCTi^PCQ (T(i„i) - Ti) + 

(48) 

U A (To - Ti ) -I- q 

Ti = F {^i) (49) 


The above symbols are analogous to those defined in Chapter II, 
Equation 48 is valid as long as the water temperature is 
below saturation, but once saturation is reached. Equation 49 
is valid. The water temperature will vary around the fluid 
circuit (see Figure 1) from an unsaturated state before the 
preheater to fully saturated and boiling within the vapori- 
zation chamber. With saturated water and vapor bubbles 
leaving the vaporization chamber and moving through the 
porous bed, heat transfer losses could cause the water temper- 
ature to drop below saturation. In order to ensure that this 
condition is detected, a saturation test is made at the begin- 
ning of each of ten sections of the porous bed. Saturation 
tests are also made in each section of tubing, within the 
accTimulator, and in the preheater and vaporization chamber. 
Equation control is needed to implement the saturation test, 
and from the results of the saturation test, employ the 
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appropriate equation in the simulation* 

As implied from the above discussion, the saturation 
test is actually two tests. One test is made to determine 
when water makes the transition from the “unsaturated to the 
saturated state” (OTSS), and a second test has to be made to 
determine when the transition is made in reverse direction, 
i.e*, from the “saturated state to the unsaturated state" 
(SSTU), 

The UTSS test is accomplished by comparing the actual 
water temperature with the saturation temperature as computed 
by a function generator which has as its input the absolute 
pressure as calculated during the simulation* The pertinent 
Fortran statements for the UTSS test are given below. These 
are followed with an explanation of the coding, 

117 RO = AFGEN(TF1S,P2S) 

122 IF(T1,GE,R0,) GO TO 8 

8 CONTINUE 

128 XU *= RO - T1 

129 T1 = M0DINT(T100,XU,XU,T1D) 

The above sequence of coding will implement the UTSS 
test. It works as follows. Assume that temperature Tl 
(Statement Number 129) has reached saturation temperature. 

The variable, RO, is the generated reference saturation 
temperature. When Tl^ RO then the program goes to 
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Statement Numl>er 8 and then proceeds to Statement Ntimher 12S, 
The control variable, KU, will have a value of zero at satu- 
ration« Statement Ntunber 129 makes use of a special CSMP 
function, MODINT, which is a “mode-controlled integrator*' 

(MCX}« The MCI functions as a regular integrator if the 
algebraic sign of the control variable, in this case XU, is 
positive. If the algebraic sign of the control variable 
becomes negative or zero, the MCI will hold the value of 
the last output until the control variable sign becomes 
positive. Thus when a saturation condition exists, the 
control variable, XU, is zero and so the MCI will hold the 
last temperature output as is desired. 

Two methods of making the SSTU test were developed to 
work with the unforced and forced temperature differential 
eguationso Equation 48 is a typical forced temperature 
differential equation. If the forcing function, q, is. zero 
then Equation 48 becomes unforced and takes the following 
forms 

P V C = P C Q(T(i_i) - Ti) + U A(To - Ti) (50) 

t 

The above symbols are analogous to those defined in Chapter II, 

The first method applies to Equation 50 and utilizes the 
knowledge that once saturation is reached the temperature 
remains constant and, therefore, the derivation of temperature 
with respect to time of Equation 50 must equal zero. Equation 
50 then becomes an algebraic equation, Equation 51, and can be 
solved explicitly for the water temperature T^, it was found 
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that Equation 51 and Equation 49 compared within three 
degrees and so Equation 51 could he used to detect a 
change from a saturated state to the unsaturated state. 

Also it app'ears that the expressions used for the heat- 
transfer coefficient of Equation 26 are correct, based 
on the results of the comparison. The pertinent Fortran 
statements which exemplify the first method of making the 
SSTU test arc given helow. This is followed with an 
explanation of the coding, 

Ti e p C Q T(i^l)/(P V C+U A) t u A To /{P V GW A) (51) 

The above symbols are analogous to those defined in Chapter II, 

117 RO = AFGEN(TF1S,P2S) 

o 

o 

122 IF<Ti.GE.RO) GO TO 8 

8 CONTINUE 

124 IF(T1 ,LT,RO )T1 = TlS 

o 

9 CONTINUE 

128 XU = RO - T1 

129 T1 = MODINT(T100,XU,XU,T1D) 

The explanation of the above sequence of coding which 
implements the first method of the SSTU test is now given. 

It is assumed that the temperature Tl has been in a 
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saturated state and enough heat transfer has taken place 
so that TlS is just dropping below RO, i.e,, temperature Tl 
is making the transition from saturation to unsaturation. 

The sequence of events begins with Statement Number 124 and 
with TlS less than RO, Since the conditional expression of 
Statement Number 124 has been satisfied, then the value of 
temperature Tl is replaced by the value of temperature TlS, 
Therefore, the control variable XU, as calculated in 
Statement Number 128, has a positive algebraic sign which 
causes the MCI of Statement Number 129 to switch from the 
hold-mode to normal integration. With the return of the 
MCI to normal integration operation, the first method of the 
SSTU test is complete. 

The development of a second method of making a SSTU 
test was motivated by the inability of the first method to 
function with equations of the form of Equation 48, When 
Equation 48 is manipulated into the form of Equation 51, it 
was found that there was an intolerable difference between 
Equation 52 and the true saturation temperature given by 
Equation 49, The symbols used in Equation 52 are similar 
to those defined in Chapter II, Equation 52 has a large 
discontinuity at the time of switching. It is highly 


Ti = P C Q T(i.x)/P V C+U A) -i- U A Tq/(P V C+U A) + 

(52) 


q/(P V C+U A) 


improbable that all the water in a particular section would 
reach saturation and instantaneously conpletely vaporize 
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into super-heated steam. This V 70 uld require an iEtipulss 
of power which is incongruent with the physical system. 

An explanation which used physical and mathematical 
arguments will toe given. The general solution to 
Equation 48 would contain the complementary or transient 
solution and the particular or steady-state solution. A 
mathematical steady-state condition is usually denoted toy 
the derivative terms vanishing and the solution tahJ.ng the 
form of the forcing function. The vanishing of the deriva- 
tive term in Equation 48, at saturation, is imposed toy the 
physical constraints of the process and not due to a steady- 
state condition in the usual mathematical sense, ThufS the 
value of the dependent variable jumped from the imposed 
steady-state value to that imposed toy the forcing function. 
The second method of making the SSTU test was made to toe 
independent of equations of the form of Equation 52. The 
pertinent Fortran coding which implements uhe secona methoa 
of performing the SSTU test is given toelovr. An explanation 


of the coding is 

then given. 

95 

REP2 


AFGEN(TF1S, PVC) 

170 

TWGD 

• 

• 

• o o 

176 

CMI2 


RBP2 - TWC 

193 

TWC 

• 

M0DINT(TWC00,CM12, CMI2, TWCD) 

194 

TWC 

• 

o 

FCNSW(CMI2, REF2, REF2, TWC) 

The explanation of ■ 

the above sequence of coding v;hich 


implements the second method of the SSTU test is now given. 
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It is assumed that the temperature TWO is saturated and 
vaporization is occurring « The vaporization and conden- 
sation rates directly affect the absolute pressure around 
the fluid circuit and the absolute pressure PVC of this 
particular section determines the saturation temperature , 
Suppose that the absolute pressure PVC were to increase „ 

Then the algebraic sign of the control-variable CMI2 would 
become positive* The MCI would resume normal integration 
of the differential equations* Statement Number 194 
utilizes the CSMP function called Function-Switch* The 
control— variable CMI2 controls the output of Statement 
Number 193* If the algebraic sign of the control-variable 
CMI2 is negative or zero, then the value of TWC is replaced 
by the value of REF2* If CMI2 has a positive sign, then the 
output of the Function-Switch is equal to TWC* Therefore, 
once saturation has been reached and vaporization is occur- 
ring, TWC cannot drop below saturation unless all heat input 
is turned off or an increase in the local absolute pressure 
occurs such that REF2 becomes greater than the temperature 
TWC, 

The UTSS test described in this chapter is used where 
appropriate throughout the simulation* The second method of 
making the SSTU test is used only with the preheater and 
vaporization sections of the system* The first method of the 
SSTU test is used in all other sections of the fluid circuit, 

A flow-diagram of the entire simulation program is given 
in Figures 13(a) through 13 (r). 
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FIGURE 13(a) o COMPUTER SIMULATION FLOW-DIAGRAM 





FIGURE 13(b), COMPUTER SliMULATION FLOW-DIAGRAM 








FIGURE 13 <c). COMPUTER SIMULATION FLOW-DIAGRAM 















FIGURE 13(e). COMPUTER SIMULATION FLOW-DIAGRAM 
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FIGURE 13(g). COMPUTER SIMUIATION PLOW-DIAGRAM 
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FIGURE 13(j)* COMPUTER SIMULATION FLOW-DIAGRAM 














FIGURE 13(k), COMPUTER SIMULATION ?*LQW-DIAGRaM 


















FIGURE 13 (m). COMPUTER SIMULATION FLOW-DIAGRAM 







FIGURE 13(n), COMPUTER SIMULATION FLOW-DIAGRAM 












FIGURE 13(p), COMPUTER SIMULATION FLOW-DIAGRAM 





FIGURE 13 (q), COMPUTER SIMULATION PLOW-DIAGRAM 
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FIGORB 13 (r). COMPUTER SIMULATION FLOW-DIAGRAM 








CHAPTER IV 


DISCUSSION AND DISP3LAY OF SIMULATION RESULTS 

The system was simulated under several operating 
conditions, i,e,, several flow- rates, a number of 
different temperature initial conditions, and various 
values of heat inputs Only tv/o runs, runs A and B are 
included. These runs are representative of the many 
runs that were simulated, i«e*, they approach the 
experiment conditions in which two-phase flow exist in 
the porous bed. The simulation runs included in this 
thesis show the response of the system for different 
temperature initial conditions throughout the system. 

The system contains 15 gallons of water. In lains 
A and B the flow-rate was 1 gpm? therefore, it will take 
15 minutes for all the water in the system to be circu- 
lated one time. In run A the system initial temperature 
conditions were lower than run B, The simulation time 
of run A was 0,8 minutes which does not allow the water 
in the fluid circuit to completely circulate. Therefore, 
run B which started at higher initial temperatures allows 
the study of system operation at a later time than run A, 
Thus run A gives the view of the system approaching the 
desired experiment conditions while run B more closely 
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approaches this condition. 

Run A 

Conditions under which run a was made are given: 

Heat input; Q1 = 10 BTU/min or zero depending upon 
the action of the preheater control. Q2 = 40 BTU/min 
or 50 BTU/min depending on the action of the vaporiza- 
tion control. 

Flow-rates started from zero flow at time equal zero, 
and achieved 1 gpm in .0025 minutes. 

Initial pressures; were equal to the absolute pressure 
in the accumulator, PAG, which was 15.834 Ibf/in^, 

Simulation times 0,8 minutes 

Initial temperatures; water throughout the system was 
212°F. 

Porous media; 208°F 
Channels 210°F 
Insulation: 75°F 

Ambients 70°F 

The results of run A are given in Figures 14 through 27, 
Under the given run conditions for run A, Figure 15 shows 
that TWG has saturated, and vaporization is occurring as 
seen in Figure 24. However, no vapor exists past the first 
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section of the porous bed since vapor cannot exist in 
those sections in which the surrounding water is not 
saturated. See Figures 18 through 21, In order for 
Vapor to exist through the porous bed the system has 
to operate a longer time so that the temperature of the 
water in the porous bed is brought up to saturation 
temperature , 

«/ 

Initially the water temperature dropped in every 
section of the porous bed. Since the initial temperature 
of each section of the porous media was set at a lower 
temperature than the surrounding water, it absorbs heat 
from the water, lowering the temperature of the water 
until heated water has flowed into the particular section. 
For example, at a flow~rate of 1 gpm, water leaving the 
vaporization chamber and proceeding through the bed would 
arrive at sections 1 through 10 ins ,0167 min,/ ,0278 
min,; ,0389 min,? ,05 min,; ,0612 min,; .0725 min,; ,0834 
min,? ,0945 min.? ,1055 min,? .1162 min., respectively. 

In 0,8 minutes none of the observable sections of 
the porous bed reached saturation temperature, a satura- 
tion condition exists if the temperature curve Intersects 
the reference saturation temperature curve. It appears 
that an intersection between the temperature curve, TF4, 
and the reference temperature curve, F2, is inttminent. See 
Figure 18, 

It is noted from Figure 25 that the pressure in the 
accumulator had reached a constant value during the 0,8 



minutes of run time. This indicates the vaporization and 

condensation have become equal? however, the heated water 

has not been circulated back to the preheater section* 

/ 

When the heated water returns to the preheater section, 
and its temperature and the temperature of the water of 
the previous cycle are not equal, then a temperature 
transient condition will exist. 

Run B 

The initial temperature conditions for run B are 
higher than run A and represent the system operation at 
a later time than run A, 

Conditions under which run B was made are givens 

Heat inputs same as run A, 

Flow~rates same as run A. 

Initial pressures; same as run A. 

Simulation times 0,9 minutes 

Initial temperatures s Water temperatures; 

TW = 214,25®F, TWC = 214.20°F, TF2 = 214,15°F, 

TF4 - 214.1QOF, TF6 = 214.08Of, TF8 = 214.05OF, 

TPIO = 214.0OF, TF12 - 213.98°F, TF14 = 213.95°F, 
TF16 213.90°F, TF18 = 213.88°F, TF20 = 213.85°F, 
TW2 = 213.83°F, TWAC = 213*80°F, TW3 = 213.780F, 

T1 s 213.77°F. 


These temperature initial conditions were set to ensure 
the temperature in any particular section of the system 
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was less than or equal to the saturation temperature for 
that section o 

Porous media s 216®F 

Channel: 216°P 

Insulation: 200*^F 

Ambient s 7 0®F 

The results of run B are given in Figures 28 through 
46, With the system initial conditions as given above, 
vaporization occurred almost immediately as noted by the 
following: (1) the intersection of the temperature curve, 

TWC, and the reference temperature curve, REF2o See 
Figure 29, (2) the increase of pressure, PAG, which is 

due to vapor displacing v/ater into the accumulator. See 
Figure 44, 

There was no initial decrease of wauer temperature, 
within the porous bed, TF2 through TP18, since the porous 
media was given a higher initxal temperature than the water. 
See Figures 32 through 40. 

During the run time, 0,9 minutes, the first four 
sections of the porous bed became saturated as the heated 
V7ater and vapor travel down the porous bed. See Figures 
32 through 35. It also appears that sections five and six 
of the porous bed are going to saturate. See Figures 36 and 
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The temperatures, TWAG and TW3, appear to have 
reached a constant value indicating the accumulator is 
receiving sufficient heat from the incoming flow to 
offset the heat transfer losses incurred in the acciamu- 
latoXo See Figures 41 and 42, 
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FIGURE 16 o RUN A, TIME RESPONSE OF VAPORIZATION HEAT-RATE 1NPUT„ 



FIGURE 17 « RUN A, TIME RESPONSE OF PREHEATER HEAT-RATE INPUT. 
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FIGURE 18 o RUN A, TIME RESPONSE OP TF4 AND P2. 
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FIGURE 20^ RUN A, TIME RESPONSE OF 


saturation reference 
temperature for section 
six of pox'ous bed 


temperature of water for 
section six of porous bed 




Degree Fahrenheit 



FIGURE 21 o RUN h, TIME RESPONSE OF TF16 AND F8„ «> 
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FIGURE 22. RUN a; TIME RESPONSE OF TWAC AND REF4 
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FIGURE 27, RUN A, TIME RESPONSE OF QUAL2, 
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RUN B, TIME RESPONSE OF TW AND REFTEM 
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FIGURE 34c RUN B, TIME RESPONSE OF TF6 AND F3 
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FIGURE 37 o RUN B, TIME RESPONSE OF TF12 AND F6 
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FIGURE 40. RUN TIME RESPONSE OF TF18 AND F9 
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FIGURE 42 o RUN B, Til® RESPONSE OF TvS3 AND RSF5 
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FIGURE 43 „ RUN B, TIME RESPONSE OF THE VAPOR CONTENT OF ADM 
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FIGURE 44. RUN B, TIME RESPONSE OF PAC 




FIGURE 45. RUN B, TIME RESPOJ^E OF QUALl 
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CHAPTER V 


CONCLUSIONS AND RECOMMENDATIONS 

An idealized mathematical model of the proposed 
flight experiment was developed in Chapter II. The 
resulting system of differential and algebraic equations 
was then prepared in a form that was representative of the 
proposed system and was digitally simulated utilizing IBM's 
Continuous System Modeling Program, The simulation program 
and several of the programming techniques were developed in 
Chapter III, In Chapter IV the foundation for later 
recommendations was laid. 

Throughout the thesis an effort was made to maintain 
a system's viewpoint while also developing the capability 
for studying the proposed system in considerable detail. 

From the simulation of this system several comments 
can be made with regard to digital simulation utilizing 
CSMP. 

CSMP is capable of handling systems which require a 
complex mathematical formulation to describe its functional 
characteristics. The implementing of physical constraints 
which are a function of the system state are relatively 
easy to incorporate into a CSMP simulation. For example, 
in this simulation it was required to swirch-in and switch- 
out the appropriate equations depending upon whether the 
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system w^s in a saturated state or an unsaturated state. 

The price pa^d for these capabilities is in total computa- 
tion time. For example, in order to obtain 0,9 minutes of 
real time simulation, it required 50 minutes of machine 
computation time. The mathematical description used to 
describe the system utilized 51 integrators and 58 function 
generators. It is well known that the digital computer is 
notoriously slow in performing numerical integration. The 
user of CSMP has six different integration methods from 
which to select. These integration methods differ in degree 
of conplexity, and consequently, in the degree of accuracy 
obtainable, and the time required to perform the numerical 
integraWon, Therefore, if computation time is of prime 
importance, the mathematical description should be simpli- 
fied, i.e,, reduce the nxamber of integrators in the simu- 
lation, Also, the use of the simpler integration methods 
should be investigated. The final selection of integration 
method will most likely depend upon a trade off between 
accuracy and computation time. 

In the system simulated, the channel and the accumu- 
lator were the only portions of the system not thermally 
insulated. Due to the requirements of the experiment, the 
channel cannot be insulated, but no such restriction is 
placed on the acc\amulator. Therefore, the heat losses 
sustained in the simulated system may be considered to be 


maximum 
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Based upon all the simulation runs and upon calculations 
made using 60 BTU/min heat input rat© (maximum available), 
the proposed system is capable of producing and maintaining 
vapor in the porous bed, 

A check was made of the coefficients of the temperature 
differential equations. For example, at saturation the 
derivative term on the left hand side of Equation 48 is 
constrained to be equal to zero. Then Equation 48 becomes 
an algebraic equation and can be solved explicitly for the 
temperature of the water in terras of the equation coeffi~ 
cients and other system temperatures. Also at saturati.on 
the water temperature is a function of pressure, and is 
given by Equation 49, The results obtained from comparing 
Equations 48 and 49 under saturated conditions compared 
with 3®F. This would indicate the coefficients of Equation 
48 and similar equations are essentially correct. 

From the system equations as given in Chapter II, and 
the results of the simulation runs, the following comments 
and recommendations are givens 

1, The system should be preheated under pressure on 
the ground and sealed, because of the 60 BTU/min 
maximiim heat input rate. By heating the system 
under pressure more heat energy can be stored. 

This is done to compensate for heat losses occur- 
ring between the time the experimental package is 
sealed aboard the space craft and time for the 
first experiment to be run. 



2, The system is extremely sensitive to pressure 
changes due to vaporization. It may be necessary 
to include a separate pressure control system. 

This feature may possibly be incorporated into 
the vaporization control system after redefining 
the control law for this control system, 

3, The design of the accumulator for this system is 
an important consideration. A trade off must be 
made between a highly pressure sensitive system 
and the values obtained in the measurement of the 
average system quality, and the use of pressure 
feedback to the vaporization control. If the 
accumulator is properly designed, then the method 
of vaporization control proposed in this thesis is 
feasible, 

4, If possible, several flow-rates in the range of 

,5 — 1,0 gpm should be selected for the experiment. 
This would allow the designer of the preheater 
control system to set the dead-band around the 
thermodynamic equilibriimi poinus as determined 
from the system simulation. 

The digital simulation of the idealized mathematical 
model is sufficiently general to allow the system designers 
to subject the simulated system to various operating condi- 
tions as they see desirable. This would cerc.ainly aid them 
in the selection of physical coiiponents and the design of 
the various control sysuems. 
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